// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_GENERAL_MATRIX_VECTOR_H
#define EIGEN_GENERAL_MATRIX_VECTOR_H

namespace Eigen {

namespace internal {

enum GEMVPacketSizeType
{
	GEMVPacketFull = 0,
	GEMVPacketHalf,
	GEMVPacketQuarter
};

template<int N, typename T1, typename T2, typename T3>
struct gemv_packet_cond
{
	typedef T3 type;
};

template<typename T1, typename T2, typename T3>
struct gemv_packet_cond<GEMVPacketFull, T1, T2, T3>
{
	typedef T1 type;
};

template<typename T1, typename T2, typename T3>
struct gemv_packet_cond<GEMVPacketHalf, T1, T2, T3>
{
	typedef T2 type;
};

template<typename LhsScalar, typename RhsScalar, int _PacketSize = GEMVPacketFull>
class gemv_traits
{
	typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;

#define PACKET_DECL_COND_PREFIX(prefix, name, packet_size)                                                             \
	typedef                                                                                                            \
		typename gemv_packet_cond<packet_size,                                                                         \
								  typename packet_traits<name##Scalar>::type,                                          \
								  typename packet_traits<name##Scalar>::half,                                          \
								  typename unpacket_traits<typename packet_traits<name##Scalar>::half>::half>::type    \
			prefix##name##Packet

	PACKET_DECL_COND_PREFIX(_, Lhs, _PacketSize);
	PACKET_DECL_COND_PREFIX(_, Rhs, _PacketSize);
	PACKET_DECL_COND_PREFIX(_, Res, _PacketSize);
#undef PACKET_DECL_COND_PREFIX

  public:
	enum
	{
		Vectorizable = unpacket_traits<_LhsPacket>::vectorizable && unpacket_traits<_RhsPacket>::vectorizable &&
					   int(unpacket_traits<_LhsPacket>::size) == int(unpacket_traits<_RhsPacket>::size),
		LhsPacketSize = Vectorizable ? unpacket_traits<_LhsPacket>::size : 1,
		RhsPacketSize = Vectorizable ? unpacket_traits<_RhsPacket>::size : 1,
		ResPacketSize = Vectorizable ? unpacket_traits<_ResPacket>::size : 1
	};

	typedef typename conditional<Vectorizable, _LhsPacket, LhsScalar>::type LhsPacket;
	typedef typename conditional<Vectorizable, _RhsPacket, RhsScalar>::type RhsPacket;
	typedef typename conditional<Vectorizable, _ResPacket, ResScalar>::type ResPacket;
};

/* Optimized col-major matrix * vector product:
 * This algorithm processes the matrix per vertical panels,
 * which are then processed horizontaly per chunck of 8*PacketSize x 1 vertical segments.
 *
 * Mixing type logic: C += alpha * A * B
 *  |  A  |  B  |alpha| comments
 *  |real |cplx |cplx | no vectorization
 *  |real |cplx |real | alpha is converted to a cplx when calling the run function, no vectorization
 *  |cplx |real |cplx | invalid, the caller has to do tmp: = A * B; C += alpha*tmp
 *  |cplx |real |real | optimal case, vectorization possible via real-cplx mul
 *
 * The same reasoning apply for the transposed case.
 */
template<typename Index,
		 typename LhsScalar,
		 typename LhsMapper,
		 bool ConjugateLhs,
		 typename RhsScalar,
		 typename RhsMapper,
		 bool ConjugateRhs,
		 int Version>
struct general_matrix_vector_product<Index,
									 LhsScalar,
									 LhsMapper,
									 ColMajor,
									 ConjugateLhs,
									 RhsScalar,
									 RhsMapper,
									 ConjugateRhs,
									 Version>
{
	typedef gemv_traits<LhsScalar, RhsScalar> Traits;
	typedef gemv_traits<LhsScalar, RhsScalar, GEMVPacketHalf> HalfTraits;
	typedef gemv_traits<LhsScalar, RhsScalar, GEMVPacketQuarter> QuarterTraits;

	typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;

	typedef typename Traits::LhsPacket LhsPacket;
	typedef typename Traits::RhsPacket RhsPacket;
	typedef typename Traits::ResPacket ResPacket;

	typedef typename HalfTraits::LhsPacket LhsPacketHalf;
	typedef typename HalfTraits::RhsPacket RhsPacketHalf;
	typedef typename HalfTraits::ResPacket ResPacketHalf;

	typedef typename QuarterTraits::LhsPacket LhsPacketQuarter;
	typedef typename QuarterTraits::RhsPacket RhsPacketQuarter;
	typedef typename QuarterTraits::ResPacket ResPacketQuarter;

	EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void run(Index rows,
														Index cols,
														const LhsMapper& lhs,
														const RhsMapper& rhs,
														ResScalar* res,
														Index resIncr,
														RhsScalar alpha);
};

template<typename Index,
		 typename LhsScalar,
		 typename LhsMapper,
		 bool ConjugateLhs,
		 typename RhsScalar,
		 typename RhsMapper,
		 bool ConjugateRhs,
		 int Version>
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE void
general_matrix_vector_product<Index,
							  LhsScalar,
							  LhsMapper,
							  ColMajor,
							  ConjugateLhs,
							  RhsScalar,
							  RhsMapper,
							  ConjugateRhs,
							  Version>::run(Index rows,
											Index cols,
											const LhsMapper& alhs,
											const RhsMapper& rhs,
											ResScalar* res,
											Index resIncr,
											RhsScalar alpha)
{
	EIGEN_UNUSED_VARIABLE(resIncr);
	eigen_internal_assert(resIncr == 1);

	// The following copy tells the compiler that lhs's attributes are not modified outside this function
	// This helps GCC to generate propoer code.
	LhsMapper lhs(alhs);

	conj_helper<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs> cj;
	conj_helper<LhsPacket, RhsPacket, ConjugateLhs, ConjugateRhs> pcj;
	conj_helper<LhsPacketHalf, RhsPacketHalf, ConjugateLhs, ConjugateRhs> pcj_half;
	conj_helper<LhsPacketQuarter, RhsPacketQuarter, ConjugateLhs, ConjugateRhs> pcj_quarter;

	const Index lhsStride = lhs.stride();
	// TODO: for padded aligned inputs, we could enable aligned reads
	enum
	{
		LhsAlignment = Unaligned,
		ResPacketSize = Traits::ResPacketSize,
		ResPacketSizeHalf = HalfTraits::ResPacketSize,
		ResPacketSizeQuarter = QuarterTraits::ResPacketSize,
		LhsPacketSize = Traits::LhsPacketSize,
		HasHalf = (int)ResPacketSizeHalf < (int)ResPacketSize,
		HasQuarter = (int)ResPacketSizeQuarter < (int)ResPacketSizeHalf
	};

	const Index n8 = rows - 8 * ResPacketSize + 1;
	const Index n4 = rows - 4 * ResPacketSize + 1;
	const Index n3 = rows - 3 * ResPacketSize + 1;
	const Index n2 = rows - 2 * ResPacketSize + 1;
	const Index n1 = rows - 1 * ResPacketSize + 1;
	const Index n_half = rows - 1 * ResPacketSizeHalf + 1;
	const Index n_quarter = rows - 1 * ResPacketSizeQuarter + 1;

	// TODO: improve the following heuristic:
	const Index block_cols = cols < 128 ? cols : (lhsStride * sizeof(LhsScalar) < 32000 ? 16 : 4);
	ResPacket palpha = pset1<ResPacket>(alpha);
	ResPacketHalf palpha_half = pset1<ResPacketHalf>(alpha);
	ResPacketQuarter palpha_quarter = pset1<ResPacketQuarter>(alpha);

	for (Index j2 = 0; j2 < cols; j2 += block_cols) {
		Index jend = numext::mini(j2 + block_cols, cols);
		Index i = 0;
		for (; i < n8; i += ResPacketSize * 8) {
			ResPacket c0 = pset1<ResPacket>(ResScalar(0)), c1 = pset1<ResPacket>(ResScalar(0)),
					  c2 = pset1<ResPacket>(ResScalar(0)), c3 = pset1<ResPacket>(ResScalar(0)),
					  c4 = pset1<ResPacket>(ResScalar(0)), c5 = pset1<ResPacket>(ResScalar(0)),
					  c6 = pset1<ResPacket>(ResScalar(0)), c7 = pset1<ResPacket>(ResScalar(0));

			for (Index j = j2; j < jend; j += 1) {
				RhsPacket b0 = pset1<RhsPacket>(rhs(j, 0));
				c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 0, j), b0, c0);
				c1 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 1, j), b0, c1);
				c2 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 2, j), b0, c2);
				c3 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 3, j), b0, c3);
				c4 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 4, j), b0, c4);
				c5 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 5, j), b0, c5);
				c6 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 6, j), b0, c6);
				c7 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 7, j), b0, c7);
			}
			pstoreu(res + i + ResPacketSize * 0, pmadd(c0, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 0)));
			pstoreu(res + i + ResPacketSize * 1, pmadd(c1, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 1)));
			pstoreu(res + i + ResPacketSize * 2, pmadd(c2, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 2)));
			pstoreu(res + i + ResPacketSize * 3, pmadd(c3, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 3)));
			pstoreu(res + i + ResPacketSize * 4, pmadd(c4, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 4)));
			pstoreu(res + i + ResPacketSize * 5, pmadd(c5, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 5)));
			pstoreu(res + i + ResPacketSize * 6, pmadd(c6, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 6)));
			pstoreu(res + i + ResPacketSize * 7, pmadd(c7, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 7)));
		}
		if (i < n4) {
			ResPacket c0 = pset1<ResPacket>(ResScalar(0)), c1 = pset1<ResPacket>(ResScalar(0)),
					  c2 = pset1<ResPacket>(ResScalar(0)), c3 = pset1<ResPacket>(ResScalar(0));

			for (Index j = j2; j < jend; j += 1) {
				RhsPacket b0 = pset1<RhsPacket>(rhs(j, 0));
				c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 0, j), b0, c0);
				c1 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 1, j), b0, c1);
				c2 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 2, j), b0, c2);
				c3 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 3, j), b0, c3);
			}
			pstoreu(res + i + ResPacketSize * 0, pmadd(c0, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 0)));
			pstoreu(res + i + ResPacketSize * 1, pmadd(c1, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 1)));
			pstoreu(res + i + ResPacketSize * 2, pmadd(c2, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 2)));
			pstoreu(res + i + ResPacketSize * 3, pmadd(c3, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 3)));

			i += ResPacketSize * 4;
		}
		if (i < n3) {
			ResPacket c0 = pset1<ResPacket>(ResScalar(0)), c1 = pset1<ResPacket>(ResScalar(0)),
					  c2 = pset1<ResPacket>(ResScalar(0));

			for (Index j = j2; j < jend; j += 1) {
				RhsPacket b0 = pset1<RhsPacket>(rhs(j, 0));
				c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 0, j), b0, c0);
				c1 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 1, j), b0, c1);
				c2 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 2, j), b0, c2);
			}
			pstoreu(res + i + ResPacketSize * 0, pmadd(c0, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 0)));
			pstoreu(res + i + ResPacketSize * 1, pmadd(c1, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 1)));
			pstoreu(res + i + ResPacketSize * 2, pmadd(c2, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 2)));

			i += ResPacketSize * 3;
		}
		if (i < n2) {
			ResPacket c0 = pset1<ResPacket>(ResScalar(0)), c1 = pset1<ResPacket>(ResScalar(0));

			for (Index j = j2; j < jend; j += 1) {
				RhsPacket b0 = pset1<RhsPacket>(rhs(j, 0));
				c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 0, j), b0, c0);
				c1 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + LhsPacketSize * 1, j), b0, c1);
			}
			pstoreu(res + i + ResPacketSize * 0, pmadd(c0, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 0)));
			pstoreu(res + i + ResPacketSize * 1, pmadd(c1, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 1)));
			i += ResPacketSize * 2;
		}
		if (i < n1) {
			ResPacket c0 = pset1<ResPacket>(ResScalar(0));
			for (Index j = j2; j < jend; j += 1) {
				RhsPacket b0 = pset1<RhsPacket>(rhs(j, 0));
				c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 0, j), b0, c0);
			}
			pstoreu(res + i + ResPacketSize * 0, pmadd(c0, palpha, ploadu<ResPacket>(res + i + ResPacketSize * 0)));
			i += ResPacketSize;
		}
		if (HasHalf && i < n_half) {
			ResPacketHalf c0 = pset1<ResPacketHalf>(ResScalar(0));
			for (Index j = j2; j < jend; j += 1) {
				RhsPacketHalf b0 = pset1<RhsPacketHalf>(rhs(j, 0));
				c0 = pcj_half.pmadd(lhs.template load<LhsPacketHalf, LhsAlignment>(i + 0, j), b0, c0);
			}
			pstoreu(res + i + ResPacketSizeHalf * 0,
					pmadd(c0, palpha_half, ploadu<ResPacketHalf>(res + i + ResPacketSizeHalf * 0)));
			i += ResPacketSizeHalf;
		}
		if (HasQuarter && i < n_quarter) {
			ResPacketQuarter c0 = pset1<ResPacketQuarter>(ResScalar(0));
			for (Index j = j2; j < jend; j += 1) {
				RhsPacketQuarter b0 = pset1<RhsPacketQuarter>(rhs(j, 0));
				c0 = pcj_quarter.pmadd(lhs.template load<LhsPacketQuarter, LhsAlignment>(i + 0, j), b0, c0);
			}
			pstoreu(res + i + ResPacketSizeQuarter * 0,
					pmadd(c0, palpha_quarter, ploadu<ResPacketQuarter>(res + i + ResPacketSizeQuarter * 0)));
			i += ResPacketSizeQuarter;
		}
		for (; i < rows; ++i) {
			ResScalar c0(0);
			for (Index j = j2; j < jend; j += 1)
				c0 += cj.pmul(lhs(i, j), rhs(j, 0));
			res[i] += alpha * c0;
		}
	}
}

/* Optimized row-major matrix * vector product:
 * This algorithm processes 4 rows at once that allows to both reduce
 * the number of load/stores of the result by a factor 4 and to reduce
 * the instruction dependency. Moreover, we know that all bands have the
 * same alignment pattern.
 *
 * Mixing type logic:
 *  - alpha is always a complex (or converted to a complex)
 *  - no vectorization
 */
template<typename Index,
		 typename LhsScalar,
		 typename LhsMapper,
		 bool ConjugateLhs,
		 typename RhsScalar,
		 typename RhsMapper,
		 bool ConjugateRhs,
		 int Version>
struct general_matrix_vector_product<Index,
									 LhsScalar,
									 LhsMapper,
									 RowMajor,
									 ConjugateLhs,
									 RhsScalar,
									 RhsMapper,
									 ConjugateRhs,
									 Version>
{
	typedef gemv_traits<LhsScalar, RhsScalar> Traits;
	typedef gemv_traits<LhsScalar, RhsScalar, GEMVPacketHalf> HalfTraits;
	typedef gemv_traits<LhsScalar, RhsScalar, GEMVPacketQuarter> QuarterTraits;

	typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;

	typedef typename Traits::LhsPacket LhsPacket;
	typedef typename Traits::RhsPacket RhsPacket;
	typedef typename Traits::ResPacket ResPacket;

	typedef typename HalfTraits::LhsPacket LhsPacketHalf;
	typedef typename HalfTraits::RhsPacket RhsPacketHalf;
	typedef typename HalfTraits::ResPacket ResPacketHalf;

	typedef typename QuarterTraits::LhsPacket LhsPacketQuarter;
	typedef typename QuarterTraits::RhsPacket RhsPacketQuarter;
	typedef typename QuarterTraits::ResPacket ResPacketQuarter;

	EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void run(Index rows,
														Index cols,
														const LhsMapper& lhs,
														const RhsMapper& rhs,
														ResScalar* res,
														Index resIncr,
														ResScalar alpha);
};

template<typename Index,
		 typename LhsScalar,
		 typename LhsMapper,
		 bool ConjugateLhs,
		 typename RhsScalar,
		 typename RhsMapper,
		 bool ConjugateRhs,
		 int Version>
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE void
general_matrix_vector_product<Index,
							  LhsScalar,
							  LhsMapper,
							  RowMajor,
							  ConjugateLhs,
							  RhsScalar,
							  RhsMapper,
							  ConjugateRhs,
							  Version>::run(Index rows,
											Index cols,
											const LhsMapper& alhs,
											const RhsMapper& rhs,
											ResScalar* res,
											Index resIncr,
											ResScalar alpha)
{
	// The following copy tells the compiler that lhs's attributes are not modified outside this function
	// This helps GCC to generate propoer code.
	LhsMapper lhs(alhs);

	eigen_internal_assert(rhs.stride() == 1);
	conj_helper<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs> cj;
	conj_helper<LhsPacket, RhsPacket, ConjugateLhs, ConjugateRhs> pcj;
	conj_helper<LhsPacketHalf, RhsPacketHalf, ConjugateLhs, ConjugateRhs> pcj_half;
	conj_helper<LhsPacketQuarter, RhsPacketQuarter, ConjugateLhs, ConjugateRhs> pcj_quarter;

	// TODO: fine tune the following heuristic. The rationale is that if the matrix is very large,
	//       processing 8 rows at once might be counter productive wrt cache.
	const Index n8 = lhs.stride() * sizeof(LhsScalar) > 32000 ? 0 : rows - 7;
	const Index n4 = rows - 3;
	const Index n2 = rows - 1;

	// TODO: for padded aligned inputs, we could enable aligned reads
	enum
	{
		LhsAlignment = Unaligned,
		ResPacketSize = Traits::ResPacketSize,
		ResPacketSizeHalf = HalfTraits::ResPacketSize,
		ResPacketSizeQuarter = QuarterTraits::ResPacketSize,
		LhsPacketSize = Traits::LhsPacketSize,
		LhsPacketSizeHalf = HalfTraits::LhsPacketSize,
		LhsPacketSizeQuarter = QuarterTraits::LhsPacketSize,
		HasHalf = (int)ResPacketSizeHalf < (int)ResPacketSize,
		HasQuarter = (int)ResPacketSizeQuarter < (int)ResPacketSizeHalf
	};

	Index i = 0;
	for (; i < n8; i += 8) {
		ResPacket c0 = pset1<ResPacket>(ResScalar(0)), c1 = pset1<ResPacket>(ResScalar(0)),
				  c2 = pset1<ResPacket>(ResScalar(0)), c3 = pset1<ResPacket>(ResScalar(0)),
				  c4 = pset1<ResPacket>(ResScalar(0)), c5 = pset1<ResPacket>(ResScalar(0)),
				  c6 = pset1<ResPacket>(ResScalar(0)), c7 = pset1<ResPacket>(ResScalar(0));

		Index j = 0;
		for (; j + LhsPacketSize <= cols; j += LhsPacketSize) {
			RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j, 0);

			c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 0, j), b0, c0);
			c1 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 1, j), b0, c1);
			c2 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 2, j), b0, c2);
			c3 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 3, j), b0, c3);
			c4 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 4, j), b0, c4);
			c5 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 5, j), b0, c5);
			c6 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 6, j), b0, c6);
			c7 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 7, j), b0, c7);
		}
		ResScalar cc0 = predux(c0);
		ResScalar cc1 = predux(c1);
		ResScalar cc2 = predux(c2);
		ResScalar cc3 = predux(c3);
		ResScalar cc4 = predux(c4);
		ResScalar cc5 = predux(c5);
		ResScalar cc6 = predux(c6);
		ResScalar cc7 = predux(c7);
		for (; j < cols; ++j) {
			RhsScalar b0 = rhs(j, 0);

			cc0 += cj.pmul(lhs(i + 0, j), b0);
			cc1 += cj.pmul(lhs(i + 1, j), b0);
			cc2 += cj.pmul(lhs(i + 2, j), b0);
			cc3 += cj.pmul(lhs(i + 3, j), b0);
			cc4 += cj.pmul(lhs(i + 4, j), b0);
			cc5 += cj.pmul(lhs(i + 5, j), b0);
			cc6 += cj.pmul(lhs(i + 6, j), b0);
			cc7 += cj.pmul(lhs(i + 7, j), b0);
		}
		res[(i + 0) * resIncr] += alpha * cc0;
		res[(i + 1) * resIncr] += alpha * cc1;
		res[(i + 2) * resIncr] += alpha * cc2;
		res[(i + 3) * resIncr] += alpha * cc3;
		res[(i + 4) * resIncr] += alpha * cc4;
		res[(i + 5) * resIncr] += alpha * cc5;
		res[(i + 6) * resIncr] += alpha * cc6;
		res[(i + 7) * resIncr] += alpha * cc7;
	}
	for (; i < n4; i += 4) {
		ResPacket c0 = pset1<ResPacket>(ResScalar(0)), c1 = pset1<ResPacket>(ResScalar(0)),
				  c2 = pset1<ResPacket>(ResScalar(0)), c3 = pset1<ResPacket>(ResScalar(0));

		Index j = 0;
		for (; j + LhsPacketSize <= cols; j += LhsPacketSize) {
			RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j, 0);

			c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 0, j), b0, c0);
			c1 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 1, j), b0, c1);
			c2 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 2, j), b0, c2);
			c3 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 3, j), b0, c3);
		}
		ResScalar cc0 = predux(c0);
		ResScalar cc1 = predux(c1);
		ResScalar cc2 = predux(c2);
		ResScalar cc3 = predux(c3);
		for (; j < cols; ++j) {
			RhsScalar b0 = rhs(j, 0);

			cc0 += cj.pmul(lhs(i + 0, j), b0);
			cc1 += cj.pmul(lhs(i + 1, j), b0);
			cc2 += cj.pmul(lhs(i + 2, j), b0);
			cc3 += cj.pmul(lhs(i + 3, j), b0);
		}
		res[(i + 0) * resIncr] += alpha * cc0;
		res[(i + 1) * resIncr] += alpha * cc1;
		res[(i + 2) * resIncr] += alpha * cc2;
		res[(i + 3) * resIncr] += alpha * cc3;
	}
	for (; i < n2; i += 2) {
		ResPacket c0 = pset1<ResPacket>(ResScalar(0)), c1 = pset1<ResPacket>(ResScalar(0));

		Index j = 0;
		for (; j + LhsPacketSize <= cols; j += LhsPacketSize) {
			RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j, 0);

			c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 0, j), b0, c0);
			c1 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i + 1, j), b0, c1);
		}
		ResScalar cc0 = predux(c0);
		ResScalar cc1 = predux(c1);
		for (; j < cols; ++j) {
			RhsScalar b0 = rhs(j, 0);

			cc0 += cj.pmul(lhs(i + 0, j), b0);
			cc1 += cj.pmul(lhs(i + 1, j), b0);
		}
		res[(i + 0) * resIncr] += alpha * cc0;
		res[(i + 1) * resIncr] += alpha * cc1;
	}
	for (; i < rows; ++i) {
		ResPacket c0 = pset1<ResPacket>(ResScalar(0));
		ResPacketHalf c0_h = pset1<ResPacketHalf>(ResScalar(0));
		ResPacketQuarter c0_q = pset1<ResPacketQuarter>(ResScalar(0));
		Index j = 0;
		for (; j + LhsPacketSize <= cols; j += LhsPacketSize) {
			RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j, 0);
			c0 = pcj.pmadd(lhs.template load<LhsPacket, LhsAlignment>(i, j), b0, c0);
		}
		ResScalar cc0 = predux(c0);
		if (HasHalf) {
			for (; j + LhsPacketSizeHalf <= cols; j += LhsPacketSizeHalf) {
				RhsPacketHalf b0 = rhs.template load<RhsPacketHalf, Unaligned>(j, 0);
				c0_h = pcj_half.pmadd(lhs.template load<LhsPacketHalf, LhsAlignment>(i, j), b0, c0_h);
			}
			cc0 += predux(c0_h);
		}
		if (HasQuarter) {
			for (; j + LhsPacketSizeQuarter <= cols; j += LhsPacketSizeQuarter) {
				RhsPacketQuarter b0 = rhs.template load<RhsPacketQuarter, Unaligned>(j, 0);
				c0_q = pcj_quarter.pmadd(lhs.template load<LhsPacketQuarter, LhsAlignment>(i, j), b0, c0_q);
			}
			cc0 += predux(c0_q);
		}
		for (; j < cols; ++j) {
			cc0 += cj.pmul(lhs(i, j), rhs(j, 0));
		}
		res[i * resIncr] += alpha * cc0;
	}
}

} // end namespace internal

} // end namespace Eigen

#endif // EIGEN_GENERAL_MATRIX_VECTOR_H
